#ifndef PLANT_SPDE2D_H
#define PLANT_SPDE2D_H

#include "Signal.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

class plant_SPDE2D:public block
{
  int    mode;
  int    LatticeSize;
  double *Km2[4];

  // States
  double *z1;
  double *z2;
  double *z3;
  double *z4;

  double T;
  double W;
  double c;
  double c2;
  double sigma;
  double **lambda;

  double dt_con; // dt_con is the dt used in numerical integration. 
  // block::dt is the time interval for output.
  double sqrt_dt_con;
  double dt_conservative_factor;

  // Random number generator supported by GSL
  const gsl_rng_type * noise_type;
  gsl_rng * noise;

  double eq_c(double T,double W);
  double eq_c2(double T,double W);
  double eq_sigma2(double T,double W);

  inline int index2D(int m,int n){return m*(mode+1)+n;};

  inline double phi1mn(int m,int n,int ix,int iy){
    return 2.0*sin(2*m*pi/LatticeSize*ix)*sin(2*n*pi/LatticeSize*iy)/pi;
  };
  //double phi2mn(int m,int n,int ix,int iy);
  inline double phi2mn(int m,int n,int ix,int iy){
    double coeff;
    if((m == 0 || m*2==LatticeSize)&&(n==0 || n*2==LatticeSize)){
      coeff = 1.0/pi;
    }
    else if(m!=0 && 2*m!=LatticeSize && n!=0 && 2*n!=LatticeSize){
      coeff = 2.0/pi;
    }
    else{
      coeff = SQRT_TWO/pi;
    }
    return coeff*cos(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);
  };

  inline double phi3mn(int m,int n,int ix,int iy){
    double coeff;
    if(n==0||n*2 == LatticeSize){
      coeff = SQRT_TWO/pi;
    }else{
      coeff = 2.0/pi;
    }
    return coeff*sin(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);
  };

  inline double phi4mn(int m,int n,int ix,int iy){
    double coeff;
    if(m==0||2*m==LatticeSize){
      coeff = SQRT_TWO/pi;
    }else{
      coeff = 2.0/pi;
    }
    return coeff*cos(2*m*pi/LatticeSize*ix)*sin(2*n*pi/LatticeSize*iy);
  };
public:
  plant_SPDE2D(int M,int L,double idt);
  ~plant_SPDE2D();
  virtual void update(double sysTime);
  virtual void reset();
  void update(double sysTime,double T,double W);
  void output(string filename);
  double getM2();
  double getR2();
};

void h2mode_2D(double *h,int L,double *z,int M);
void mode2h_2D(double *z,int M,double *h,int L);
double randn();
#endif